
% This is the comparison of two methods, using the simulation.
% The function is intended to reproduce the Figure 5 in the paper.
% Finished on June 15th 2021. Based on Equation (35).

function [T, Mean, Std, ExpectedProfit, profit_pro, profit_rrall, Histo_5,  profit_pmax, best_i, ...
    best_ip, Summary, Num_in_slot_p, Percentage_in_slot_p, Figure_9, Figure_16, ImproveRange] = Simu_Figure5 (LINE)

    global P 
    P = [.25, .5, .75];

    tic

    % Import data 'LINE.mat' that contains many simulation lines.
    [Num_simu, N] = size(LINE);
    ImproveRange = zeros (1, Num_simu);
    ExpectedProfit = zeros (1, Num_simu);
    best_ip = zeros (Num_simu, N);   % Store the best i_slot for proposed ones.

    for t = 1:Num_simu
        fprintf("\t fig5 -> %d/%d \n", t, Num_simu);
        % Firstly, for the proposed policy.
        line = LINE(t,:);
        [ best_i, profit_pro, ~, ~, ~ ] = prop_policy(line, N);
        [ profit_pmax, situation ] = max(profit_pro);
        best_ip(t, 1: situation) = best_i(1:situation);        % We can get the detailed information of best_ip.
        ExpectedProfit (1, t) = profit_pmax;
        % Then, find the corresponding Round-Robin Policy.
        [~, profit_rrall,~ ,~ ,~ ] = round_robin (line, N);
        profit_rmax = profit_rrall (1, situation);
        % Calculate the improvement range!
        ImproveRange (1,t) = 100* (profit_pmax - profit_rmax) / profit_pmax;
    end

    % Then, histogram ImprovementRange!
    figure(5)
    Histo_5 = histogram(ImproveRange, 50);
    xlabel ('Percentage Difference')
    ylabel ('Frequency Distribution')
    title ('Fig 5: Improvement over Round Robin.')
    Mean = mean(ImproveRange);
    Std = std(ImproveRange);
    T = toc;

    % For Figure 9:

    best_ip( best_ip==0)=[ ];
    Summary = tabulate(best_ip(:));

    Num_in_slot_p = zeros(1, 8);
    Percentage_in_slot_p = zeros(1, 8);

    for n = 1:8
        Num_in_slot_p(1, n) = Summary (n, 2);
        Percentage_in_slot_p(1,n) = Summary (n, 3);
    end

    figure(9)
    Figure_9 = plot(Percentage_in_slot_p, 'ko-');
    xlabel ('Slot')
    ylabel ('Average % Patients Assigned')
    title ('Fig 9: Percentage Assignments per slot.')

    figure(16)
    Figure_16 = histfit (ExpectedProfit,150);
    axis([1240 1340 0 80]);
    xlabel ('Expected Profit')
    ylabel ('Frequency Distribution')
    title ('Fig 16: Frequency Histogram of Expected Profit for 2500 Sequences.')

end
